## negative log likelihood function
get.y.prob.from.MLE = function(data, theta.coef, eta.coef,gop.coef,phi.coef, cnt, rep=500, seed){
  
  # if(parn %in% names(pars)){
  #   pars[[parn]][1:length(pars[[parn]])] = par.update
  # } 
  
  L0 = cbind(rep(1,nrow(data)),data[,1:pL0])
  l0 = L0 %*% 2^(1:(pL0+1))
  L0 = unique(L0,margin = 1)
  l0.index = match(l0, unique(l0))
  Nb = nrow(L0)
  
  theta.coef = matrix(theta.coef[c(1,3:(pL0+2), 2:(pL0+2))],,2)
  colnames(theta.coef) = c("l0","l1")
  theta = exp(L0 %*% theta.coef)
  
  phi.coef   = matrix(phi.coef[c(1,4:(pL0+3),2,4:(pL0+3),
                                      3:(pL0+3))],,3)
  colnames(phi.coef)   = c("a0","a1","a.") 
  
  eta.coef        = matrix(eta.coef[c(1,6:(5+pL0), 2, 6:(5+pL0), 3, 6:(5+pL0),
                                      4, 6:(5+pL0), 5:(5+pL0))],,5)
  colnames(eta.coef) = c("00","01","10","11",".")
  phi   = exp(L0 %*% phi.coef)
  eta   = expit(L0 %*% eta.coef)
  gop   = exp(L0 %*% gop.coef)
  
  ### Step 1: get ratio between two "adjacent" probabilities
  ratio = t(sapply(1:Nb, function(i) GetRatio(theta[i,], phi[i,],eta[i,])))
  dim(ratio) = c(Nb, 5, 2)
  dimnames(ratio) = list(1:Nb, c("ak","lk","aK","lK","l0"),c("0","1"))
  
  ### Step 2: get the maximal ratio against the "baseline p000" 
  #   for each observed unit (only depends on L0)
  k.max = GetMaxKGen(ratio, K, Nb)
  k.max[k.max == Inf] = .Machine$double.xmax
  check = which(k.max < 1e10)
  
  ### Step 3: get the largest p (only depends on L0)
  p.max = rep(1 - 0.1^10, Nb)
  if(length(check)>0)     p.max[check] = PMaxEstRough(ratio[check,,], K, k.max[check], gop[check],rep = rep, seed = seed)
  
  
  ### Step 3: get the ratio against the largest p for each observed unit
  ratio = ratio[l0.index,,]
  p.max = p.max[l0.index]
  k.max = k.max[l0.index]
  p.y.fn <- function(A, L){
    # A = data[,(pL0+1):(pL0+K+1)];   colnames(A) = 0:K
    # L = data[,(pL0+K+2):(pL0+2*K + 2)];   colnames(L) = 0:K
    colnames(A) = 0:K;  colnames(L) = 0:K
    k = GetKGen(A, L, ratio, k.max, K)
    k[k.max == .Machine$double.xmax] = 0
    
    p.y   = p.max * k
    return(p.y)
  }
  
  return(p.y.fn)
}

y.logit.est <- function(data, cnt, K){
  
  N    = nrow(data)
  Y = data[,"Y"]
  Ak   = data[,(pL0+K)]
  Lk   = data[,(pL0+2*K+1)]
  AkLk00 = as.numeric(Ak == 0 & Lk == 0)
  AkLk01 = as.numeric(Ak == 0 & Lk == 1)
  AkLk10 = as.numeric(Ak == 1 & Lk == 0)
  AkLk11 = as.numeric(Ak == 1 & Lk == 1)
  
  X = data[,1:pL0]
  # temp   = list()
  # for(i in 1:(K+1)) temp[[i]] = X.base
  # X      = do.call(rbind, temp)
  
  X.model = model.matrix(~ AkLk00 + AkLk01 + AkLk10 + AkLk11  + X - 1)
  wt      = cnt
  
  a.model = glm(Y~ AkLk00 + AkLk01 + AkLk10 + AkLk11  + X - 1,
                data=data.frame(cbind(X.model, Y)), family=binomial(), weights = wt)
  
  # Finished training
  # Now want to evaluate outcome when setting Ak and Lk to user specified value
  p.y.fn <- function(A, L){
    # A = data[,(pL0+1):(pL0+K+1)];   colnames(A) = 0:K
    # L = data[,(pL0+K+2):(pL0+2*K + 2)];   colnames(L) = 0:K
    
    Ak   = A[,K+1]
    Lk   = L[,K+1]
    AkLk00 = as.numeric(Ak == 0 & Lk == 0)
    AkLk01 = as.numeric(Ak == 0 & Lk == 1)
    AkLk10 = as.numeric(Ak == 1 & Lk == 0)
    AkLk11 = as.numeric(Ak == 1 & Lk == 1)
    pred.model.matrix = model.matrix(~ AkLk00 + AkLk01 + AkLk10 + AkLk11  + X - 1)
    pred.prob <- expit(predict(a.model, newdata=data.frame(pred.model.matrix)))
    return(pred.prob)
  }
  
  return(p.y.fn)
  
}


prob.score.func <- function(b.coef,L,X){
  b.cov <- cbind(rep(1, length(L)), L, X)
  return(expit(b.cov %*% b.coef))
}


p.L.func <- function(eta.coef, ak, lk, X){
  eta.coef <- eta.coef[c(1:4,6:9)] # ignore gamma. because we don't use it
  AkLk00 = as.numeric(ak == 0 & lk == 0)
  AkLk01 = as.numeric(ak == 0 & lk == 1)
  AkLk10 = as.numeric(ak == 1 & lk == 0)
  AkLk11 = as.numeric(ak == 1 & lk == 1)
  pred.model.matrix = model.matrix(~ AkLk00 + AkLk01 + AkLk10 + AkLk11  + X - 1)
  # pred.prob <- expit(predict(a.model, newdata=data.frame(pred.model.matrix)))
  pred.prob <- expit(pred.model.matrix %*% eta.coef)
  return(pred.prob)
}
